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We critically discuss the application of the Wertheim’s theory to classes of complex associating 
fluids that can be today engineered in the laboratory as patchy colloids and to the prediction 
of their peculiar gas-liquid phase diagrams. Our systematic study, stemming from perturbative 
version of the theory, allows us to show that, even at the simplest level of approximation for 
the inter-cluster correlations, the theory is still able to provide a consistent and stable picture 
of the behavior of interesting models of self-assembling colloidal suspension. We extend the 
analysis of a few cases of patchy systems recently introduced in the literature. In particular, 
we discuss for the first time in detail the consistency of the structural description underlying 
the perturbative approach and we are able to prove a consistency relationship between the 
valence as obtained from thermodynamics and from the structure for the one-site case. A 
simple analytical expression for the structure factor is proposed. 

Keywords: Colloidal suspensions, Wertheim thermodynamic perturbation theory. 
Associating fluids. Structure of fluids. 


1. Introduction 


Recently, there have been interesting developments of techniques for the synthesis 
of new colloidal patchy particles in the laboratory [l|], including seeded growth, 
swelling, and phase separation. Whereas in the laboratory relatively less work has 
been done on the thermodynamic characterization of self-assembly of these par¬ 
ticles, from a theoretical point of view, or in recent computer experiments, these 
kind of associating fluids and their clustering and phase behavior are actively 
studied M. 

In principle, statistical mechanics should be able to describe all equilibrium 
phases. However, the strong and confined attractions responsible of association 
call for a more clever approach than brute force. In particular, it has been found 
useful to describe an associating fluid as one where there are ric species of clusters 
made of a number i of particles, denoted i-mers. Many dehnitions of cluster are 
possible [iMi either of a geometric nature or of a topological one, depending on 


the spatial arrangement of the bonded particles. If we measure the concentrations 
of the z-mers in an associating fluid we will hnd that they are functions of the ther¬ 
modynamic state: For one-component systems, the temperature T and the density 
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p of the fluid. Then, special statistical mechanics approaches have been developed 
to obtain such information and phase diagrams from models of interactions. 

In our previous work we compared two theories for cluster equilibria, the 
Wertheim association theory [iMi and the Bjerrum-Tani-Henderson theory [13- 
[ 2 ^ and we showed that for Uc = 2 the two approaches coincide when inter-cluster 
correlation are ignored, i.e. the system behaves as an ideal gas of clusters. Nonethe¬ 
less, the simple and elegant perturbation theory described in Wertheim’s work, is 
able, unlike the one of Bjerrum-Tani-Henderson, to describe the case of tt-c —> oo 
fluids. Due to this fact, Wertheim theory is able to describe the liquid phase, thus 
giving access to the study of liquid-gas coexistence in a coherent way, while the 
Bjerrum-Tani-Henderson one is not. The first order in the Wertheim perturbation 
theory approximation is a simple but very useful tool. At high temperature, the 
associating fluid reduces to the “reference” fluid that can also be considered as the 
one obtained from the associating fluid switching off all attractions. However, in 
its original form, the theory is only applicable when some “steric incompatibility” 
conditions are fulhlled by the associating fluid: A single bond per site, no more 
than one bond between any two particles, and no closed loop, or ring, of bonds. 

Patchy colloids are systems of current experimental and theoretical 0,113 inter¬ 
est. Simple models for their interactions, for example fluids of hard-spheres deco¬ 
rated with attractive sites distributed on their surface, are well suited for applica¬ 
tion of Wertheim theory. For particles with M identical bonding sites, Bianchi et 
al. 0-0] discovered the “empty liquid” scenario as M approaches two, i.e. when the 
clusters allowed in the fluid are just the “chains”. Even more rich phenomenology 
is found when there are sites of two different kinds 00 and “junctions” forma¬ 
tion becomes possible. Such structures become responsible for a re-entrance of the 
liquid branch of the binodal, and for “rings” formation 0,0 . Moreover, extending 
Wertheim theory beyond its steric incompatibility conditions, the rings formation 
has been found to be responsible for a re-entrance also in the gas branch and the 
appearance of a second lower critical point (recently appeared studies which further 


extend Wertheim theory to allow also for doubly bonded sites [28|-l30||). From all 


these studies emerged how Wertheim theory has very good semi-quantitative agree¬ 
ment with exact Monte Carlo simulations, when applied to these one-component 
patchy particle fluids (especially so at the level of the clusters concentrations behav¬ 
ior). Far from being a purely theoretical speculation, these fluids can be engineered 
in the laboratory []| from patchy colloids. 

In the present work, while critically reviewing such theoretical results, in par¬ 
ticular elucidating the role of the accuracy of inter cluster correlations, we will 
discuss the solution of the Wertheim theory applied to hard-spheres with M iden¬ 
tical bonding sites and with sites of two different kinds. Our analysis is intended to 
be as simple and systematic as possible while re-analyzing the many works found 
in the literature on various particular highly idealized associating colloidal suspen¬ 
sion models. This will allow us to treat the ring forming systems of Rovigatti et al. 
fully analytically as freeR jointed chains. We show that also the results in Ref. 

3l[ |. extending Russo et al. [a, l7(| results to take into acccount the “X-junctions” 
formation, and in particular the existence of charcteristic “R” shaped spinodals are 
largely independent on the choice of the reference system correlations. Moreover, 
we find indication of a gas-liquid coexistence with a critical point at extremely low 
densities and temperatures at r < 1/3, with r the ratio between the gain in energy 
between the bond of two unlike sites and the one between two like sites. 

We also study in detail the relationship between structural and thermodynamic 
information within Wertheim theory, and in particular between the effective valence 
as obtained from the thermodynamics and from the structure. 
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The paper is organized as follows: In Section [2] we introduce the thermodynamic 
quantities we will take under consideration in the rest of the work; in Section [3] we 
will review Wertheim association theory in the light of the present work needs, the 
problem of identical attractive site (Section I3.1.2h . and the problem of attractive 
sites of two different kinds lSection l3.1.3p : in Section [3]2] we introduce the problem 
of the gas-liquid coexistence; in Section 13.31 we comment on the relevance of the 
pair-potential microscopic level of description; in Section [4] we sistematically re¬ 
analyze many results obtained applying Wertheim theory to specific fluids with 
identical sites fSection l4.1h and sites of two different kinds (Section l4.2p . We show, 
in a systematic way, that all the results present in the literature are structurally 
stable with respect to changes in the reference system accuracy; in Section [6] we 
determine a simple analytical expression for the radial distribution function which 
we then use to calculate the valence; in Section [3 we determine a simple analytical 
expression for the structure factor; Section [8] is for final remarks. 


2. Thermodynamics 

Consider a one-component fluid of N associating hard-sphere (HS) particles in a 
volume V at an absolute temperature T = l/jSkB with ks Boltzmann constant. 

The Helmholtz free energy H of a hard-sphere associating fluid can be written 
as a sum of separate contributions 


^ -^0 T ^mf T -^bondi 


( 1 ) 


where Aq is the free energy of a hard-sphere fluid at a density p = N/V, is 
the mean-held contribution due to the dispersion forces, and is the change 

in the free energy due to association. We will generally use the notation a{p, T) = 
a = A/N for the free energy per particle. 

The hard-sphere free energy per particle in excess of the ideal gas one is accu¬ 
rately given by the Carnahan and Starling expression 




Arj — ‘irf 

TT^’ 


( 2 ) 


where p = (7r/6)/9(j^ is the packing fraction of the hard-spheres of diameter a. So 
that adding the ideal gas contribution j3aid = ln(pA^/e), with A the de Broglie 
thermal wavelength, we obtain oq = aid + Oq^. 

The mean-held contribution has the van der Waals form 


where the constant e^/ is the measure of the strength of the mean-held attractions. 
The addition of this contribution to Aq is essential to have a gas-liquid coexistence. 

From a microscopic point of view one can see, for example, the mean held con¬ 
tribution as arising from the hrst order in /3 in a high temperature expansion of a 
thermodynamic perturbation theory treatment of the square-well (SW) huid, with 
the HS taken as the reference system. So, the free energy of the corresponding 
associating huid will be given by A = Asw + Abond- But, as we will see in Section 
m one can have gas-liquid coexistence with just A = Aq + Ahond for a properly 
chosen Abond- 
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We can define a unit of length, S, and a unit of energy, £, so that we can introduce 
a reduced density, p* = pS^, and a reduced temperature, T* = ksTjE. 

The association contribution A^cmd will be discussed in the next section. 


3. Associating fluids 

We recall here the main result of Wertheim association theory [Mi. We write 
the bond free energy per particle abond such that the full free energy per particle of 
the associating fluid can be written as a = oq + abond-, where oq is the contribution 
of the reference fluid, the one obtained from the associating fluid setting to zero 
all the bonding attractions. We discuss the importance of the choice of a proper 
pair-potential for the fulfillment of the steric incompatibility conditions in the mi¬ 
croscopic description of the fluid. And we discuss the problem of the determination 
of the gas-liquid coexistence line (the binodal) in our one-component fluid. 


3.1. Wertheim statistical thermodynamic theory 

In Wertheim theory [Mi one assumes that each hard-sphere of the one- 
component fluid is decorated with a set T of M attractive sites. Under the as¬ 
sumptions of: [i.] a single bond per site, [ii.] no more than one bond between any 
two particles, and [hi.] no closed loop, or ring, of bonds, one can write in a first 
order thermodynamic perturbation theory framework, valid at reasonably high 
temperatures, 


PaZnd = Yi 


osF 


2 


M 

T’ 


( 4 ) 


where Xa = N^/N is the fraction of sites a that are not bonded. We will also 
introduce the symbol Xi to denote the concentration of clusters made of a number 
i of particles. We will always use a Greek index to denote a specific site. We can 
solve for the x^ from the “law of mass action” 


1 


Xo = 


1 + pXj/Jsr XpE^aji 


a 6 r 


( 5 ) 


where the probability to form a bond, once the available sites of the two particles 
are chosen, is given by piS.ap = pAp^ and approximated as 


Aop = J 5o(’'i2)(/o/3(12)>ni,n2dri2. 


( 6 ) 


Here go is the radial distribution function of the reference system, fap is the Mayer 
function between site a on particle 1 and site (5 on particle 2 (see Section l3..Sjl , 
and {.. denotes an angular average over all orientations of particles 1 and 

2 at a fixed relative distance ri 2 . Eq. ([5]) should be solved for the real physically 
relevant solution such that limp^o = 1- Even if we cannot exclude the possibility 
of having multiple solutions satisfying to this condition we never encountered such 
a case in the present work. Clearly we cannot assign any physical value to the 
branches with Xa ^ [0,1]. 

At high temperatures A^p —> 0 and ^ 1, which means we have complete dis¬ 
sociation. At low temperatures (Wertheim theory is a high temperature expansion 




May 27, 2015 Molecular Physics MolPhys-associationRev 


but here we just mean the formal low T limit of the first order Wertheim results) 
00 and Xa —> 0, which means that we have complete association. 

The number of attractive sites controls the physical behavior. Models with one 
site allow only dimerization. The presence of two sites permits the formation of 
chain and ring polymers. Additional sites allow formation of branched polymers 
and amorphous systems. 

3.1.1. One attractive site 


The case of a single attractive site was carefully considered in our previous 
work [2] where a comparison between the Wertheim theory and the Bjerrum-Tani- 
Henderson theory 2[H26l| was made. 


3.1.2. Identical attractive sites 


Another simple case we can consider in Wertheim theory is the one with M 
identical attractive sites of kind A (we will always use a capital letter to denote 
a site kind). Now the law of mass action for x = xa (the fraction of unbonded 
specific sites of kind A) is solved by 


2 

1 + ^/lA^~AMpK' 


( 7 ) 


with A = 

The free energy contribution due to association is now given by 

PaZnd = M[\nx-xl2) + Ml2. ( 8 ) 


In this case xi = x^. 

3.1.3. Attractive sites of two kinds 

A more complex case in Wertheim theory is the one with Ma identical attractive 
sites of kind A and Mb identical attractive sites of kind B. Now the law of mass 
action reduces to the following system of two coupled quadratic equations 

XA + MapAaax'a + MbpAabxaxb = 1 , ( 9 ) 

xb + MbpAbbXb + AIapA.abxaxb = 1 , ( 10 ) 

which admits in general a set of 4 different solutions for {xa,xb) from which it is 
necessary to single out the physically relevant one. In the event that there is no 
attraction between a site of kind A and a site of kind B then Aab = 0 and the 
system simplifies to 


2 

1 + -y/l + 4lMaPAaA 

(11) 

2 

1 + + AMbpAbb 

(12) 


In the event that there is no attraction between sites of the same kind it simplifies 
to 


XA — 2/{l + {Mb — Ma)pAab + 


^[1 + {Mb - Ma)pAabY + 4M^/9A^s} 


( 13 ) 
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and xb obtained exchanging A B in the equation above. 

The free energy contribution due to association is now given by 

~ — XaI‘^) + MAj‘2: + 

MsilnxB - xb/2) + Mb/2. (14) 

In this case xi = x^'^x^^. 

3.2. The gas-liquid coexistence 

In order to determine the gas-liquid coexistence line (the binodal) one needs to 
find the compressibility factor = f3p/p, with p the pressure, and the chemical 
potential p of the associating fluid according to the thermodynamic relations 

(15) 

/3At(/5, r) = =z + l3a. (16) 

The coexistence line is then given by the Gibbs equilibrium condition of equality 
of the pressures and chemical potentials of the two phases 

Pgz{pg,T) = piz{pi,T), (17) 

PKPg^T) = Pp{pi,T), (18) 

from which one can find the coexistence density of the gas Pg{T) and of the liquid 
pi{T) phases. 

The critical point {pc, T^) is determined by solving the following system of equa¬ 
tions 


dzp 

dp 


Pc,Tc 



= 0, 

(19) 

= 0. 

(20) 


3.2.1. The mean field case 


For the HS fluid in the presence of just a van der Waals mean field free energy 
contribution, described by Eq. ([T]) without the last association term, the thermo¬ 
dynamics is parameter free. We take the diameter of the spheres a as the unit of 
length (so that p* e [0, -\/2] with -\/2 the close-packing reduced density) and e^/ as 
the unit of energy. Solving the Gibbs equilibrium conditions of Eqs. (I17jl - (ll8j) we 
find the binodal of Fig. [Hand from Eqs. (I19|l - (l2n|) we find the critical point. 

We can see this case as describing a thermodynamic perturbation theory approx¬ 
imation for a SW fluid to first order in fi small 35|. Monte Garlo simulations of 
the SW fluid are well known to show a gas-liquid binodal with the critical point 
shifting at lower temperatures and higher densities as the width of the attractive 
well decreas es l36l . [37 1. 

Recently |38|| it was shown through numerical simulation and theoretical ap¬ 
proaches that a binodal with two maxima, implying the existence of a low-density- 
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Figure 1. Gas-liquid binodal for the HS plus the van der Waals mean field term. The circle is the critical 
point at p* 0.249129, T* ^ 0.180155, and 2^ == 0.358956 [H. 

liquid and a high-density-liquid, can arise solely from an isotropic interaction po¬ 
tential with an attractive part and with two characteristic short-range repulsive 
distances. 

We consider the binodal of Fig. [T] as “standard” in the sense that the gas branch 
Tg{p) is a monotonously increasing function of density and the liquid branch Ti{p) 
a monotonously decreasing function of density. We will see in the next section that 
using Wertheim association theory it is possible to obtain non standard binodals by 
replacing the mean-field contribution Amf with a proper association contribution 


3.3. Microscopic description: Importance of the pair-potential 

The fluid is assumed to be made of particles interacting only through a pair- 
potential (/>(12) = (^(ri, Oi, r 2 , O 2 ) where r* and 0* are the position vector of the 
center of particle i and the orientation of particle i respectively. 

To give structure to the fluid we further assume that the particles have an 
isotropic hard-core of diameter a with 

</.(12) = 4(ri2) + $(12), (21) 

where ri 2 = |ri 2 | = |r 2 — ri| is the separation between the two particles 1 and 2 
and 


Mr) 


- 1-00 r ^ cr 
0 r > a ’’ 


( 22 ) 


The anisotropic part $(12) in Wertheim theory is generally chosen as 


$(12) = 2 S (23) 

asr per 


where 


ra/3 = r2 + d/3(n2) - ri - dQ,(ni), (24) 

is the vector connecting site a on particle 1 with site /3 on particle 2. Here d^ is the 
vector from the particle center to site a with da < The site-site interactions 
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^ 0 are assumed to be purely attractive. The Mayer functions introduced in 
Section [Q are then defined as fa^{ 12 ) = exp[—/3'0o/3(?"a^)] ~ 1- 
Wertheim theory depends on the specific form of the site-site potential only 
through the quantity Aa/s of Eq. ([6]), as long as the three conditions of a single 
bond per site, no more than one bond between any two particles, and no closed 
loop of bonds, are satisfied. A common choice, for example, is a square-well form 




€a/3 ^ ^ da/3 
0 r >dai 3 ' 


(25) 


where > 0 are site-site energy scales, the wells depths, and are the wells 
widths. In this case we must have d^ + d /3 > a — da /3 moreover we will have 

AaP = Kaf3{(T, da/3, “ !)• (26) 

We will also call limp_>o Ka /3 = Ka /3 some purely geometric factors. Remember 
that limp^offo(?’) = ©(’^ — f) with 0 the Heaviside step function. 

Another common choice is the Kern-Prenkel patch-patch pair-potential model 

3 ^. 


4. Structural stability of Wertheim theory 

There has recently been some relevant progress on the study of several complex as¬ 
sociating fluids through Monte Carlo (MC) simulations and theoretically through 
the Wertheim theory outlined above. The comparison between the two approaches 
shows semi-quantitative agreement, between the exact MC results and the approx¬ 
imated theoretical results, at the level of description of clusters concentrations and 
of gas-liquid binodal. We will here return on some of the ^sterns studied from 
Bianchi et al. [3-[i], Russo et al. [1,0, and Rovigatti et al. 0,0 from a unified per¬ 
spective, and concentrating ourselves on the structural stability of the Wertheim 
theory, i.e. we will show that all the qualitative non standard features of the phase 
diagrams at a large extent do not depend on the accuracy of description of the 
reference system. 


4.1. Identical sites 


The case of hard-spheres with a number M of identical attractive sites in various 
geometries on the surface of the spherical particle has been studied by Bianchi et al. 
[3-[5t]. They showed that the properties of the resulting fluid are largely independent 
from the sites geometry 0] . And the gas-liquid binodal has a liquid branch moving 
at lower densities as M decreases. In particular the binodal vanishes for M ^ 2, 
a scenario that they called “empty liquid”: The critical temperature Tc{M) and 
critical density pdM) are such that limyf^2 Tc = Tc > ^ and limyf^2 Pc = 0. There 
is then the formation of a homogeneous disordered material at small densities below 
Tc, i.e. a stable equilibrium gel. Moreover, in their fluid with M = 2, Bianchi et al. 
observed linear “chains” formation: “chaining”. 

This is quite different from what happens in fluids of Kern and Frenkel patchy 


hard-spheres varying the patches surface coverages [40!]. In Ref. [40( a study of 


criticality similar to the one of Bianchi was made varying the attractive patch 


surface coverage y. As the surface coverage y vanishes, lim,j,_ 
was found in such cases. 


I Tc = lim,^^o Pc = 0 
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Liu et al. [3^ repeated Bianchi study for a system of square-wells (SW), instead 
of HSs as in the Bianchi case, with a number M of identical attractive sites. In 
their study the gas-liquid coexistence remains also for M —> 0, as expected in view 
of the comments of Section 13.2.11 


4-1-I. Gas-liquid binodal 

With M identical sites of kind A we have in the site-site interaction eaa = e 
which we take as unit of energy and again we take a as unit of length. 

We now choose a = uq + cLbond with the association part given by the Wertheim 
theory Eq. ([4]) with M identical sites (see Section r3.1.2p . 

Following Ref. [3] we choose the identical sites distributed on the surface of the 
spherical particle and 


dAA = d 



all ^ 0.120U, 


(27) 


which guarantees that each site is engaged at most in one bond. Moreover we 
approximate the radial distribution function of the reference system with its zero 
density limit taking IAaa = A = — l] and using, in Eq. (I26p . the following 

expressions 


</aa(12)> 


mAA{r) 


(e^^ - l)mAA{ri2) n2 > cr, 



r)‘^{2d 

6rcj2 


a + r) , 

- a < r < a d 

r > a + d 


J 'tT+d 

mAA{r')r‘^dr 

a 

= 7r(i^(15(T -I- 4:d)/30a‘^ 
^ 0.332 X 


(28) 

(29) 


(30) 


In Fig. [2] we show the evolution of the gas-liquid binodal as a function of M, the 
only free parameter in Wertheim thermodynamic perturbation theory. Compared 
with Fig. 4 of Bianchi et al. we see how the qualitative behavior stays the same 
even if the two figures differ slightly quantitatively due to our further approxima¬ 
tion of taking the radial distribution of the reference system equal to one in the 
range where bonding occurs. This shows how the Wertheim theory is robust in its 
qualitative phase diagram predictions. The binodal appears to be always a standard 
one. And, as we can see from the figure, upon approaching M ^ 2 the coexistence 
disappears. Bianchi et al. called this phenomenon the empty liquid scenario. It 
in particular tells us that the fluid with M = 2, with the two sites chosen at the 
spherical particle poles in order to avoid the formations of rings (closed loops of 
bonds), is made only by chains and does not admit a gas-liquid coexistence. The 
non-integer M cases can be realized through a binary mixture UMM- 
From the point of view of Wertheim theory the reason for this scenario can be 
explained simply by looking at the low temperature limit for the bond contribution 
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Figure 2. (color online) Top panel: Evolution of the gas-liquid binodal as a function of M. The continuous 
thick black line is the locus of the critical points for M g] 2, 5]. Bottom panel: Pressure-temperature 
diagram. 


to the pressure 


HpZm - 

2M2Ap2 M 

=- K —>- p. ul 

(1 + VI + 4MAp) 2 

From which immediately follows that for M > 2 the pressure as a function of 
density on a low temperature isotherm shows a van der Waals loop at low densities, 
which implies the occurrence of a gas-liquid coexistence region. 


4.2. Sites of two kinds 


Tavares et al. l4^ studied the case of HS with three sites, two identical A 
sites at the poles and a third B one. In addition to chaining, here they observe 
the formation of “junctions”: “branching”; rings formation is inhibited in these 
cases since the A sites at the poles have very small well widths and the B site 
position is chosen so as to avoid small bond loops, i.e. triangular and square ar¬ 
rangements of bonded particles. Two types of junctions are possible in models 
where A A bonds are responsible for the chaining: X-shaped junctions, due to BB 
bonds, and Y-shaped junctions, due to AB bonds. They found that when two 
of the three interaction strengths vanish simultaneously, there can be no liquid- 
vapor coexistence. These correspond to the limits of non interacting linear chains 
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{(^AA ^ 0, ess = eAB = 0), dimers {esB ^ 0,6^^ = cab = 0), and hyperbranched 
polymers (cab ^ 0, eAA = ess = 0) of Eq- (fl^ . They also showed that the phase 
transition always disappears as eaa 0. Moreover they showed that whereas “X- 
junctions” only yield a critical point if their formation is energetically favorable, 
fluids with “Y-junctions” will exhibit a critical point, even if forming them raises 
the energy, provided this increase is below a certain threshold. 

Russo et al. 0,13 extended Tavares study to the case of two identical small A 
sites at the poles and nine equispaced identical big B sites on the equator. Killing 
the interaction between two B sites {cbb = 0) they observed the formation of 
chains and Y-junctions (and possibly hyperbranched polymers for ^abI^AA large 
enough) and eventually a re-entrant behavior of the liquid branch of the gas-liquid 
binodal pinched at low temperatures. 

Rovigatti et al. [1, @] extended Russo model selecting an off-pole position of the 
A sites, thus adding the possibility of “rings” formation, and observed re-entrance 
both in the gas and in the liquid branch of the binodal with a second lower critical 
point where the coexistence curves closes itself at low temperatures without the 
pinch. They needed to relax assumption [hi.] in Wertheim theory 45l-l47l|. 


4 . 2 . 1 . Gas-liquid binodal 

Russo et al. 0] studied the case of sites of two different kinds when the site-site 
interaction is restricted to cbb = 0 (no X-junctions). Then choosing as unit of 
energy caa and again a as the unit of length the Wertheim theory depends on only 
five parameters: r = cabI^aa > 0 and Ma^ Mb, KaAt Kab- 

We now choose a = oq -f abond with the association part given by the Wertheim 
theory Eq. ([4]) with sites of two different kinds (see Section I3.1.3P . In particular 
with the condition eBB = 0, Eqs. (j ^ - (fT0]) admit just a set of 3 different solutions 
for {xa,xb) from which it is necessary to single out the real physically relevant 
one such that lim^o^^A = hmp^oa^B = 1- 

Eollowing Ref. we choose Ma = 2, Mb = 9 (see Fig. [3|) and = 1.80 x 
= 1.56 X 10“^cr^. In order to fulfill the Wertheim condition [i.], of 
a single bond per site, the small A sites are meant to reside at the particle poles 
and the big B sites equispaced on the particle equator. The choice of 
and the large Mb make branching entropically favorable. We then approximate 
Aaa = - 1) and Aab = - 1). 

In Fig. 0] we show the evolution of the gas-liquid binodal as a function of r. Once 
again, comparing with Fig. 3 of Russo et al. we observe a complete qualitative 
agreement, even if in our calculation we further approximated the radial distribu¬ 
tion of the reference system equal to one independently of density. We see that 
for r < 1/2 we have a non standard binodal with a re-entrant liquid branch and a 
“pinched” shape evidence that indeed the topological phase separation of Tlusty 
and Safran [i^ is observed. Russo et al. @| were able to provide a qualitative ex¬ 
planation for this behavior by analyzing the energetic of the junction formation 
process: since the energy cost of forming a chain end is echain = eAA/2 > 0 and the 
energy cost of forming a Y-junction is ey- junction = -^AB + eAA/2 = eAA(l/2 - r), 
for r < 1/2 we have ey-juncUon < 0, and at low temperatures only chains, which 
we already saw that do not phase separate, are present. 

They are also able to conclude that phase separation occurs only if r > 1/3. For 
r < 1/3, the energy cost of forming junctions being too high or, alternatively, the 
entropy gain being too small to offset the loss of translational entropy of chains in 
the liquid phase. 

This behavior can be understood by looking at f{T,p;r) = djSp/dp = dl3{po + 
P^nd)/^P- Differently from Bianchi et al. case now we have limp^o dfdp'^.^Jdp = 0. 
The zeroes of / are two lines in the (p, T) plane, one for the minima of the pressure 
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A 



Figure 3. (color online) Pictorial view of a colloidal particle with attractive sites of two different kinds: 
Two A sites on the poles and nine B sites on the equator. 


and one for the maxima. The union of the two lines is called the spinodal line for 
the coexistence. The equal area construction tells us that the binodal line encloses 
the spinodal line and the two lines are tangent at the critical point. In Fig. [5] we 
show a tridimensional plot of / for r = 0.36, 2/5,1/2 as a function of temperature 
and density. Clearly the three different scenarios do not depend on the specihc 
values of Kaa, Kab, Ma, Mb which only influence the region in the phase diagram 
(yO, T) where we have the van der Waals loop. 

The cluster populations for the chain ends, 2xa, and Y-junctions, 9(1 — xb), 
along the binodal were studied in Ref. @] and are shown in their Fig. 4. From Fig. 
9 of Ref. we see how the mean value of the number of bonds per particle (the 
valence), 2(1 — xa) + 9(1 — xb), tends to 2 at low temperatures, i.e. the fluid tends 
to be formed essentially by chains which, in agreement with Bianchi et al. analysis, 
are unable to sustain the gas-liquid coexistence. 

The study of Russo et al. differs substantially from the Janus fluid case 22-23, 
where it is found a re-entrant gas branch for the gas-liquid binodal. 

Rovigatti et al. extended Russo study to take account of rings formation. In 
this case the expression for the Wertheim bond free energy per particle of Eq. (|14p 
with Ma = 2 should be corrected as follows 


l^oZnd = In (vXb’') - xa - ^xb + 

Mb Go 
P 


1 + 


(32) 


where Gn is the nth moment of the rings size distribution 


Gr, = 


2 i^Wii2pAAAy)\ 


(33) 


here imin is the minimum ring size, y is the fraction of particles with the two A sites 
unbonded, and IF/ is the number of configurations of a ring of size i. Assuming for 
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T* 



Figure 4. (color online) Top panel: Evolution of the gas-liquid binodal as a function of r. The continuous 
thick black line is the locus of the critical points for r g] 1/3, 1/2]. Middle panel: Pressure-temperature 
diagram. Bottom Panel: binodals of Russo et al. 0 Fig. 4 as obtained from their analysis (lines) of the 
Wertheim theory and from their MC simulations (points); the big circles are their predicted critical points. 


the rings the freely jointed chain level of description we can approximate 


{i + 


iji - 1) y (-IP f i-l-2j V-^ 
Stt 3 )'\ 2 ) 


(34) 


for I the smallest integer which satisfies I ^ (i — l)/2 — 1. Expression (j34j) is due 
to Treloar and is the value of the end-to-end distribution function for a freely 
jointed chain of i links, when the end links are the length of one link apart (the 
link length is equal to the diameter of a sphere which we take to be our unit of 
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1-0.36 



P 

1 = 1/2 



Figure 5. (color online) Tridimensional plots of /(T, p;r) = df3p/dp (green surface) for r = 0.36,2/5,1/2 
from top to bottom. Also shown is the plane f = 0 (blue surface). For r = 1/3 the two surfaces become 
tangent at small temperatures and small densities. For r > 1/2 the minimum in the pressure moves at 
larger densities at smaller temperatures. 


length). For i » 1 it has the following asymptotic behavior 

/ Q \ 3/2 

(i + l)Wi+i ^ j e-3/2^ i » 1 , (35) 

The laws of mass action of Eqs. (©-(UnD, for Ebb = 0, should now be corrected 
to take into account of the 7 ^ 0 as follows 

x\ = y{l-Gi/p), (36) 

l-XA = MbpAabxbxa + 2pAaax‘a + Gi/p, (37) 

1 - XB = 2pAabxaxb- (38) 


Note that solving for xa Eq. (1361) and for xb Eq. (I38|) and substituting into Eq. 
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()37l) one finds an equation in y only, which always admits just one solution y with 
the properties 0 ^ y ^ 1 and limT^o 2 / = 0. 

In Fig. [6] we show our theoretical numerical results for the gas-liquid binodal of 
the ring forming fluid. A comparison with Fig. 1 of Rovigatti et al. [St] shows again 
a good qualitative agreement between the two calculations. In our calculation we 
retained the first 50 terms in the convergent series of Eq. (j.S3j) and chose Mb = 9 
and Aaa, Aab as before. As we can see the rings formation is responsible for the re¬ 
entrance in both the gas and liquid branches of the binodal and for the appearance 
of a second lower critical point. At r = 0.37 we could not find a coexistence line, 
leaving a system for which self-assembly is the only mechanism for aggregation. 




0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 
T* 


Figure 6. (color online) Top panel: Evolution of the gas-liquid binodal as a function of r. The thin lines 
are the binodals of Fig. [4| The thick lines are the results obtained for the rings forming fluid. Bottom 
panel: Pressure-temperature diagram. 


In particular upon approaching the upper critical point, at T = r“, if we make 
a reversible transformation going from the liquid phase to the vapor phase on an 
isotherm, at T < T“, we will have, as usual 


AS = 




(39) 


with AS the change in entropy S = —{dA/dT)jqy, 5Q the infinitesimal heat ex¬ 
changes along the path of the transformation, the “latent” heat of vaporization, 
and m the mass of the fluid. Whereas Rovigatti et al. show that upon approach¬ 
ing the lower critical point, at T = T^, in the same transformation at T > T^, one 
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finds 


I 


T 


A,,m 


= A5 < 0, 


(40) 


so that the “latent” heat of vaporization changes sign as T varies from T“ to T^. 
This can be seen directly from our pressure-temperature diagram of Fig. [6] using 
the Clapeyron-Clausius formula 511. 

Rovigatti analysis neglects the rings with AB bonds. We think that their inclu¬ 
sions may have dramatic effects on the phase diagram. 

4-2.2. A possible extension 

It is possible to extend Russo et ai. 0,0 results allowing for the ebb AO 
condition, responsible for the X-junctions formation (h| . The analysis for just 
three sites, two of kind A and one of kind B, can be found in Refs. [431 . |44 | were, 
interestingly enough, it is found the disappearance of criticality as caa —*■ 0. In 
our extension we can introduce an additional parameter s = eBBl^AA > 0. One 
immediately verifies that the law of mass action of Eq. (©-(doD admits now 4 
solutions {xa,xb) from which one has to determine the physical one such that 
xa,xb g [0,1] and \\m.p^QXA = Omp^QXB = 1- Clearly in the limit r ^ 0 the 
problem is similar to the one of Bianchi et ai. 0] (compare Eqs. (|lip - (jl2p and 
Eq. dZI) and in the limit s —> 0 we fall back to Russo et al. 0,0 case. We are 
interested in the non-trivial case: /S.aa A Abb or Ma A Mb- We will choose 
for Ma, Mb, and the same values of the Russo’s case of Section r4.2.1l 
Moreover we will choose = K\j^. Again one has limp_»o = 0. 

For s small we are still able to see the re-entrant liquid scenario contrary to the 
predictions of Ref. 4^. In other words we are able to observe a re-entrant liquid 
branch even in the presence of X-junctions in the fluid, as long as the energy cost 
for their formation, ex-junction = ^aa{). — s), is positive and big enough. This is 
shown in Fig. [3 The figure also shows how an “R” shaped spinodal is possible in 
these cases with a majority of Y-junctions in correspondence of the coexistence 
region at high temperature, a majority of X-junctions in correspondence of the 
coexistence region at low temperature, and a majority of chains in between in 
corres pon dence of the bottleneck in the “R”, in agreement with the study of Tavares 
et ai. [^. Moreover we find gas-liquid coexistence also for r < 1/3 as long as s is 
large enough. This is shown in Fig. [8] from which it is also apparent the existence 
of a gas-liquid coexistence with a critical point at extremely low densities and 
temperatures, unpredicted by the study of Tavares et al. j3l|. As a matter of fact 
the critical temperature can be made small at will by a proper choice of the control 
parameters s; the spinodal being essentially independent from r. 


5. Break-down of the theory 

Apart from the necessity to fulfill the steric incompatibility conditions the 
Wertheim theory will break-down in the following cases: 


5.1. Low temperature limit 

Both the Wertheim theory and the canonical Monte Carlo simulation break-down 
at low temperatures. The Wertheim theory is a high temperature perturbation 
theory. The first order version that we have been using until now clearly breaks- 
down at low temperature when from the mass action law ([5|) follows that Xq ^ 0 
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Figure 7. (color online) Tridimensional plots of f{T,p]r,s) = djSp/dp (green surface) and of the plane 
f = 0 (blue surface) for (r, s) = (2/5,1/5). We show two plots one at high temperature and one at low 
temperature because the (x^,a:g) physical solution determination changes in the two regions of the phase 
diagram. The negative / in the high temperature and high density corner of the lowest plot is due to 
another change in the physical solution determination. 


which in turn produces an undefined bond free energy Q. Also the usual Monte 
Carlo simulation will break-down at very low temperatures. In fact, imagine we 
have to break a bond with a single particle move. Then the total energy difference 
between the final configuration and the initial one would be e and we would need 
around single particle moves. So at low temperatures we would need a very 

long simulation in order to fully explore configuration space. Depending from the 
computational resources at one disposal the range of inaccessible temperatures, 
before the solidification at zero temperature where the fluid chooses spontaneously 
the minimum potential energy configuration, may vary. Even if it is possible that 
patchy fluids, with short-ranged and tunable pair-interactions and with limited 
valence, will not crystallize at zero temperature remaining a liquid in that 
limit. 


5.2. Infinite number of attractive sites 

The Wertheim theory will not be applicable anymore to particles decorated with 
too many attractive sites. In the limit of an infinite number of sites uniformly 
distributed over the particle surface one recovers the square-well fluid or the mean 
field solution of Section 13.2.11 

















May 27, 2015 Molecular Physics 


MolPhys-associationRev 


1 = 1/4, s= 1/4 
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Figure 8. (color online) Tridimensional plots of f{T,p;r,s) = d^p/dp (green surface) for (r,s) = 
(1/4,1/4), (1/4,1/2), (1/8,1/2). Also shown is the plane / = 0 (blue surface). As we can see the spinodals 
of the two cases (n, s) = (1/4,1/2), (1/8,1/2) look essentially the same. 


6. The radial distribution function 


Using the fact that the angular average of the functional derivative of the free 
energy per particle respect to the angle dependent pair-potential is equal to p/2V 
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times the radial distribution function of colloid centers, we can write 

(41) 

(42) 

where we denote with {...) the orientational average, and in the second equality 
we used Eq. dH and Eq. ([21]). 

To make some progress we use the following property 

= -maf3{r) - (faf}} (43) 

where in the last equality we used Eqs. ()28|) and ()29p . Erom Eq. ([6|) follows 

6Aafs/6(fal3{l2)} = 47rri25o(?’i2)^a/3(n2), (44) 

where la/si'f') is equal to one on the support of (/q,/?) and zero otherwise. Next we 
observe that 






P \(i<^(l,2). 


. ^ 2 1 

= 9o{r) + - 


p dvrr^ 




7£r 


5x^ 




,/3er 


5Xr. 


^ [lja,/3sr 


— y 


6 x 








'M2 ^ 
o,/3£r 




dx-, 


BA 


0.(3 


(45) 


where M is the total number of sites per particle and in the last equality we used 
the chain rule. So we obtain 


g{r) = g^ir) 


1 + 


1 

APp 


E 

o,/3,7£r 


2 

1 - 

X.y 


dx^ 




(46) 


where the terms can be determined from the law of mass action, Eq. ([5|). In 
particular, using the symmetry A^p = Apa, it follows 


i y (^1 _ 


'3C(y_X 


(47) 


Erom Eq. (|46p we can extract the contact value for the radial distribution function 
g{(^^) = go{cF^) X 


E 


M2p 


Q;,/3,7Gr 


Bx^ 

x^ J BA^p 


m. 




(48) 
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where map{a) is the product of the two solid angle fractions for the bond when 
two particles are located at relative center-to-center distance a. For example for 
the Kern and Frenkel pair-potential [s^ we would have map = XaXp with Xpatch 
the patch surface coverage. In the Bianchi et al. case of Section [Q we have 
instead maa{o') = id/cr)^/3, from Eq. For go{a'^) we can use the analytic 

solution to the Percus-Yevick approximation for the hard-sphere fluid [s^, namely 

go{a+) = + (49) 


Next we observe that, since pg{r)A7rr‘^dr gives the number of particles in the spher¬ 
ical shell [r, r + dr] around a particle fixed on the origin, the coordination number 
can be estimated as follows 


■I 


r-\-d 


Cn = p\ A-nr^gQ{r) x 


1 + 


M'^p 




a,/3,7Gr 


Xr, 


dx^ 


dA, 




-map{r)ed^°‘^ 


dr, 


(50) 


where d = min{dQ,^}. The mean number of bonds per particle (the valence), vt = 
— Xq), can be also estimated from the structure as follows 


vs = Cn- lim Cn 
r^oo 


1 

1^ 


E 

a,/3,7Gr 


i-A 


dx^ 

d^aP 


^ap • 


Then using Eq. (I47p we immediately find 



Xa^P^aP 

a,peT 


m2 


asr 


(51) 


(52) 


where the last equality follows from the law of mass action, Eq. ([5]). The sought 
for consistency between the valence calculated from the thermodynamics and the 
valence calculated from the structure only holds in the single site per particle case, 

M = 1. 

Eor example, for M identical sites we find vt = M{1 — x) and, choosing Kern- 
Frenkel patches for which d represents the width of the attractive square well of 
each patch and x the patch surface coverage, from Eq. (jlT)) follows 


a, = 


r(J+d 

J ' 


^'Kr‘^go{r) 


1 -h 


dr. 


(53) 


7. The structure factor 

We then determined the structure factor S{k) = 1-1- ph{k) with h{r) = g{r) — 1 
the total correlation function and the hat denotes the Fourier transform. 
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7.1. Identical sites 

For the case of Bianchi et al. of Section O we find 

1 + x^m(r)e^^ 

sin(A:r) 

- - -r dr, 


S{k) = l + 47rp\ 


U X 


(54) 


where x is given by Eq. ([7]) and m(r) is given by Eq. (|29p . Choosing for go{r) = 
0(r — cj) the one obtained from the zero density limit of the hard-sphere fluid, we 
find the “triangular” approximation result of Eq. (lAip of Appendix Erom this 
result follows immediately 


lim S{k) = 1-1- 

k^O 


20r/ - 8Mr?(e^" - 1)) (15d^ + Ad^) - 

4 ( 5 -h VEy/b + 4d‘^M'q{eP^ - 1)(15 -h Ad) 


5 -h -h Ad^Mriie^^ - 1)(15 -h 4d)^ . 


Moreover we find 


lim S'(O) = 1 - 8r] + —, 
T^o ^ ^ ' M 


^im S'(O) = 1 — 8r] + ( -I- -m ) ij, 


(55) 

(56) 

(57) 


whereas for the structure factor of the reference system we have S'o(O) = 1 — 8r/. 
In Fig. [9] we show the structure factor of Eq. ()A1I) for M = A and T* = 0.1, rj = 

0 . 1 . 



Figure 9. (color online) Structure factor for M = 4 and T* = 0.1 ,77 = 0.1 in the Bianchi et ai. case using 
for the radial distribution function of the reference system, qq, the zero density limit of the hard-sphere 
fluid. Also shown, for comparison, is the structure factor of the reference system, Soik) = 1-1- 24r)(k cos(fc) — 
sm{k))/k^. 

A comparison with the simulation results of Sciortino et al. 0] (see their Fig. 
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13) at M = 2 and T* = 0.055 shows that approximation (1551) breaks-down at high 
densities. This is shown in Fig. [10] where the data of Sciortino et al. simulations 
are compared with the isothermal compressibility sum rule, 


S(0) 


5 

dp \ dp J \ 


(58) 


and the relationship between the activity A and the density is obtained 

through Eq. m- We think that the fact that the structure as determined by 
the Eq. (j54p does not satisfy the isothermal compressibility sum rule of Eq. (|58p 
is a thermodynamical inconsistency not universally recognized for the Wertheim 
theory. In order to find accurate results for the structure one needs to solve the 
Wertheim Ornstein-Zernike equation with an appropriate closure (s^ . 



Figure 10. (color online) Structure factor at zero wave-number as a function of density for M = 2 and 
T* = 0.055 in the Sciortino et al. simulations of Ref. 01, from the thermodynamic route (T) of the 
isothermal compressibility of Eq. ilSSl . from the structure route (S) of Eq. (I55II . and from the zero wave- 
number limit of Eq. Isit taking as a reference system the Percus-Yevick analytic solution for hard-spheres 
(S PY). 


8. Conclusions 

We have critically analyzed some recent applications of the Wertheim perturba¬ 
tion theory to classes of associating fluids of with non standard phase diawams 
and increasing complexity which can be today engineered in the laboratory []|. In 
particular, we have illustrated the strong structural stability of the theory, which 
allows to get a first correct qualitative understanding of the resulting phase dia¬ 
grams, even at the simplest level where all correlations of the reference system are 
neglected. 

For fluids of hard-spheres with M identical bonding sites Bianchi et al. M 
discovered the “empty liquid” scenario as M approaches two, i.e. in the presence 
of “chains” only. The phenomenology when there are sites of two different kinds is 
more rich [^, and one can have “junctions”, responsible for a re-entrance of the 
liquid branch of the binodal, and “rings” [1, [^, responsible for a re-entrance also 
in the gas branch and the appearance of a second lower critical point. 

In our detailed analysis of these results, we show that all the important conclu¬ 
sions on the qualitative behaviour of the phase diagrams can be derived uniquely 
from theoretical analytical considerations without the need of inputs from sim¬ 
ulation results. For example, for the case of rings forming fluids we used as the 
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partition function of an isolated ring the Treloar analytic expression for a freely 
jointed chain, unlike Rovigatti et al. who use a fit of the MC data. This 

approximation makes immediately available a useful tool of analysis of complex 
phase diagrams even in absence of more accurate but heavy numerical results. 

Also in the case of the more demanding condition of the presence of X-j unctions 
we find that, when the energy gain for an X-junction formation, s, is low enough, 
we still observe a re-entrant liquid branch for r < 1/2 in the fluid, eventually with 
an “R” shaped spinodal in agreement with the study of Tavares et al. [3l| . When s 
is sufficiently large we observe gas-liquid coexistence also at r < 1/3 in agreement 
with the predictions of Ref. 4j] . In these latter cases a gas-liquid coexistence with 
a critical point at an extremely low density and temperature, unpredicted by the 
work of Tavares et al. [3lJ, can be observed. 

Moreover, we have discussed in detail the consistency between structural and 
thermodynamic description within Wertheim perturbation theory and in particular 
the valence as obtained from the thermodynamics and from the structure. We can 
conclude that while the overall structural information underlying the first order 
perturbative level is not accurate, the theory provides a consistency condition on 
the estimate of bonded particles, which is satisfied only in the one-site case. An 
analytical expression for the radial distribution function and the structure factor 
has also been proposed. 


Appendix A. The structure factor in the “triangular” approximation 

Choosing go{r) = 0(r — cr) in Eq. (l5H) with m{r) defined as in Eq. (|2^ . we find 


S{k) = 


1 -I- 80r/ 



dOd'^k^Mr] — 2Ad^k^Mri) cos{k) + 


[md'^k^Mr] + 2M^k^Mr] + I0d^k^)e^^ cos{k) + 

cos{k) + 

(— 15/c^ -I- 90d‘^k^Mri + 2Ad^k‘^Mr]) sin(fc) -I- 
{—90d'^k‘^Mr] — 2Ad^k'^Mrj)e^^ sm{k) + 

{15d‘^k‘^ + 30)6^*^ sin(A;) -I- 

—3\/5k‘^'\J 5 -I- 4d^Mri{ed^ — 1)(15 -I- Ad) sm{k) + 


30{dk cos(A:(l -I- d)) 


sin(A:(l -I- d)))e^^ 


/ 


5 + V5^5 + Ad^Mr]{eP^ - 1)(15 + Ad) 


(Al) 


From this expression one immediately sees that the high temperature limit, /3 ^ 0, 
of the structure factor is independent from the number of sites, M. 
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